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This paper introduces a one-dimensional NPZD-model developed to simulate biological 
activity in a turbulent ocean water column. The model consists of a system of coupled 
semilinear parabolic equations. An initial-boundary value problem is formulated and the 
existence of a unique positive weak solution to it is proved. The existence result is derived 
using a variational formulation, an approximate model and a fixed-point method. It is 
shown that the qualitative analysis performed still applies if different parameterizations 
of several biological processes found in the biogeochemical modeling literature are used. 
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1. Introduction 

Within the scope of global climate studies, authors carrying out modeling research 
use marine ecosystem models with increasing degrees of complexity. Complexity in 
such models can arise from the number of biological compartments, or state vari- 
ables, which are taken into account, as well as from the parameterizations used to 
model interactions between these compartments. The number of variables can vary 
from one to more than ten. At least two variables, nutrients (N) and phytoplankton 
(P) are necessary to model primary production, that is to say the transformation of 
mineral nutrients into primitive biotic material using external energy, provided by 
the sun (Taylor et at ^*). However, in order to study the ocean carbon cycle, the 
main biological processes which have to be understood and estimated are primary 
production, but also the export of organic matter from the surface to deep ocean 
layers and organic matter remineralization. The simpliest model able to represent 
all these processes contains four variables, nutrients (N), phytoplankton (P), zoo- 
plankton (Z) and detritus (D). This type of model is termed the NPZD-model and 
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different variants of it are used in numerous studies. All these models are similar 
from a structural point of view but authors use different parameterizations to model 
fluxes between biological compartments. 

Complexity can also arise from the spatial resolution of the physical dynamics to 
which biological variables arc submitted. Many model set-ups arc zcro-dimensional 
and biological variables correspond to ocean mixed-layer values (e.g., Fasham et al. 

Steele and Henderson Spitz et al. Fennel et al. ^^). Others are one- 
dimensional, considering that the ocean, in some particular places, can be modeled 
with a good approximation by a turbulent water column (e.g., Prunet et al. 
Doney et al. ®, Levy et al. Memery et al. ^^). Finally, in some studies, 
the biological model is integrated in a three-dimensional circulation model (e.g., 
Fasham et al. , Moisan et al. Levy et al. , Carmillet et al. The 
question which has motivated this work is: are all these models well-posed? Of 
course, it seems difficult to study all of them and in this work we concentrate flrst 
on a one-dimensional NPZD-model and then discuss the possible generalization of 
our result. The three-dimensional version of the model we consider is proposed 
in Levy et al. and a one-dimensional version of a similar model, containing 
six biological variables, is used by Faugeras et al. to assimilate data from the 
JGOFS-DYFAMED time-series station in the North- Western Mediterranean Sea. 

Mathematically, the biological model under consideration is a system of coupled 
parabolic semilinear equations to which initial and boundary conditions are added. 
Under certain hypotheses, this general type of initial-boundary value problem can 
be transformed to an abstract Cauchy problem and studied using the theory of 
semigroups (following Chapter 6 of the book by Fazy for example). In their paper 
Boushaba et al. * used results on semigroups to provide a mathematical analysis of 
a model describing the evolution of a single variable phytoplankton. Although the 
model they considered is three-dimensional the biological reaction terms are quite 
simple since only production and mortality of phytoplanckton are represented. The 
model we propose here seems to be more realistic and has already been numerically 
validated using observations from the DYFAMED time-series station (Levy et al. 
^"^ , Faugeras et al. ^^). 

If data are regular enough the semigroup method can enable the existence of 
classical solutions to be proved. However it does not enable parabolic equations 
with time-dependent irregular coefficients to be easily handled. Since this is the 
case in the NPZD-model we consider a variational formulation approach is more 
attractive. The main purpose of this paper is to address the issue of the existence 
of weak solutions to this particular one-dimensional model. The method we propose 
is inspired from the work of Artola ^ in which an existence result for a semilinear 
parabolic system is derived using a fixed-point argument. We introduce an approx- 
imate model and prove the existence of weak solutions to this model using this 
method. Wc then pass to the limit in the approximate model to prove the existence 
of weak solutions to the NPZD-model. Furthermore, as the variables of the model 
represent concentrations they should be positive. We show this is the case. 
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Wc shall now briefly outline the contents of the paper. In the next section we 
introduce the equations of the one-dimensional NPZD-model and give some com- 
ments on the different parameterizations used. In Section 3 we set the mathematical 
framework and state our main result, which is proved in Sections 4 and 5. The goal 
of Section 6 is twofold. First wc show that the existence and positivity results still 
hold when different parameterizations found in the literature are used. Second, we 
address the issue of uniqueness of solutions. In order to prove uniqueness we need 
the nonlinear reaction terms to satisfy a local Lipschitz condition. We show this is 
the case in our particular model. 

2. Presentation of the one-dimensional NPZD-model 
2.1. Equations of the model 

In this section wc give the equations of the one-dimensional NPZD-model and for- 
mulate the initial-boundary value problem which will be studied. 

Let us first of all justify the use of a one-dimensional model. We have in mind 
numerical studies (Faugeras et al. Levy et al. Memery et al. conducted 
with such onc-dimcnsional models. In these studies simulations arc forced with 
physical data (wind stress, heat fluxes, evaporation-precipitation) and validated 
by comparison with biogeochemical data (chlorophyll and nitrate) collected at the 
DYFAMED station. This station is located in the Northwestern Mediterranean 
Sea and is an interesting test case for several reasons. First, several biogeochemical 
production regimes that take place in the world ocean are found here. Secondly, the 
station is far enough away from the Ligurian Current to be sufficiently protected 
from lateral transport, thereby permitting a one-dimensional study. 
In the above cited numerical studies the biogeochemical model is integrated in a 
one-dimensional physical model, which simulates the time evolution of velocity, 
temperature, salinity and turbulent kinetic energy (TKE). Advection is neglected 
even though this might result in a crude approximation in summer during strong 
wind events (Andersen and Prieur ^). The only dynamic process which is taken 
into account is vertical diffusion. 

The onc-dimcnsional NPZD-model consists of four coupled semilinear parabolic 
equations. Before introducing them let us give some notations. In all the following 
we denote the nutrient, phytoplankton, zooplankton and detritus concentration 
vector by, 

C = {N, P, Z, D) = (Ci, C2, C3, Ci), 
and the reaction terms by, 

f = (/jv, /p, fz, Id) = (/i, /2, /s, fi)- 
The equations of the NPZD-model read as follows. For i = 1 to 4: 
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ot Ox Ox ox 

fir' fir' 

^(t,0) = ^(t,L) = 0, te]0,T], 

[ Ci{0,x)=CUx), xe]0,L[, 

(2.1) 

with 

' fN{t,x,C) = {-tip{l - j)Li{t,x,P)LNP + t^zZ + HdD)l]o,i]{x) 
+iT{P + Z + D))l^i^^{x), 
fp{t,x,C) = ifXp{l-^)Liit,x,P)LNP-GpZ-mpP)l^o,i]{x) 

+(-rP)l],,i[(x), 
fz{t,x,C) = {apGpZ + adGoZ - nizZ - iJzZ)lp^i]{x) 
+ i-rZ)^,L^ix), 

fD{t,x,C) = {{1 - ap)GpZ - adGoZ + mpP + mzZ - i2dD)l-^o^q{x) 

+{-rD) l^i^Llix). 

(2.2) 

T is a fixed time. In immerical simulations, system (2.1) is intregated over a period 
of time which can vary from one month to a few years. 

L is the depth of the water column under consideration {L w 1000 m), I is the 
maximum depth of the euphotic layer {I w 200 m). 

l]o,z] andl]j_i[ are the usual indicator functions, 

^ . . r 1 if a;e]0,/], 
^1m(^) = | otherwise. 

Si^4 is the Kronecker symbol, 

§ _ r 1 if i = 4, 
''^ ' otherwise. 

Ncuman boundary conditions at a; = and x = L express the fact that there is no 
flux through the surface of the ocean and through the ocean floor. 

Initial concentrations, Cf, satisfy Cf{x) > for all x g]0,L[. 



The different parameters which appear in the reaction terms are strictly pos- 
itive constants. All of them are shown in Table 1. A schematic representation of 
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the model is shown on Figure 1. Let us note that parameters 7,ap and aa satisfy 
l-7>0, l-ap>0 and 1 - > 0. 

The nonUnear functions Lj,Ln,Gp and Go are given exphcitfy in the following 
subsection, and more details about the model can be found in Levy et al. 



Table 1: Parameter values 



parameter 


name 


value 


unit 


half-saturation constant 


kn 


0.5 


mmolNm~'^ 


maximal grazing rate 


9z 


0.75 


day~^ 


half-saturation constant for grazing 




1 


mmolN.m~^ 


assimilated fraction of phytoplankton 




0.7 




assimilated fraction of detritus 


ad 


0.5 




zooplankton excretion rate 




0.1 


day~^ 


phjrtoplankton mortality rate 


rup 


0.03 


day~^ 


zooplankton mortality rate 




0.03 


day~^ 


detritus remineralization rate 


fid 


0.09 


day~^ 


detritus sedimentation speed 


Vd 


5 


m.day~^ 


maximal growth rate 




2 


day~^ 


exsudation fraction 


7 


0.05 




remineralization rate 


T 


0.05 


day~^ 



Figure 1: Schematic representation of the compartments and processes of the NPZD 
surface layer model. 
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2.2. Comments and hypotheses 

1. The mixing or diffusion coefficient, d{t,x), is obtained diagnostically from 
TKE (Caspar et al. ^*). In modeling studies it is considered in the first 
approximation that biological variables do not influence physical variables. As 
a consequence biological tracers are vertically mixed with the same coefficient 
as temperature and salinity. This coefficient is an output of the physical model 
and data for the biological model. Consequently it does not depend on C. We 
are thus given once and for all a mixing coefficient d{t, x). It strongly varies in 
space and time and we can not assume it is particularly regular (Lewandosky 
^^). The usual basic assumption made in the mathematical literature as well 
as in numerical studies is the following. We suppose that 

< rfo < d{t,x) < doo, a.e. in ]0,T[x]0,L[. 

2. Because of functions Ijo.i] and the equations of the model are not the 
same above and below the depth I which physically corresponds to the depth 
at which the action of light on the system becomes negligeable. This corre- 
sponds to a discontinuity of the reaction terms fi{t, x, C) at the point x = I. 
It is the choice of modclization made by Levy et al. . Above the depth I 
the reaction terms correspond to the schematic representation of the model 
shown on Figure 1. The basic biogeochemical fluxes are represented using a 
minimum number of prognostic variables. Nutrients allow the estimation of 
production to be made. Zooplankton mortality and dctrital sedimentation 
feed the particle export flux. Below the depth / remineralization processes 
are preponderent and the surface model does not apply. Instead decay of 
phytoplankton, zooplankton and detritus in nutrients parameterize reminer- 
alization. More details about the modeled biogeochemical processes can be 
found in Levy et al. ^"^ . In the following points we give the analytical expres- 
sion of the nonlinear terms which are used. 

3. Z/jv, Gp and Gd are nonlinear functions. 

• Ltv parameterizes the nutrient limitation on phytoplankton growth. It 

N 

follows the Michaelis-Menten kinetic, Ln — . The possible nulli- 

kn+N 

N 

fication of the term fc„ + N , invites us to define, L^r = ■ — -. This 

kn + \N\ 

formulation will be used in the following. We will show that if initial con- 
centrations are positive then concentrations always stay positive, thus the 
two formulations are equivalent. Let us remark that, 

L;v is defined and continuous on IR, 

\Ln{n)\ < 1, vAf em.. 
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• Gp and Go are the zooplankton grazing rates on phytoplankton and 
detritus. The formulation used is a squared Michaelis-Menten response 
function: 

In the remainder of this paper we use the following properties: 
Gp and Gd are defined and continuous on IR, 
|Gp| <5., VPGIR, 
\GD\<gz,yDGTR. 

4. The limitation of phytoplankton growth by light is parameterized by, 

Li{t, x,P) = l- e^p{-PAR{t, X, P)/kpar), 

kpar is a positive constant. The photosynthetic available radiation, PAR, is 
predicted from surface irradiancc and phytoplankton pigment content by a 
light absorption model according to Levy et al. Prom a biological point of 
view, the fact that PAR depends on P is important. This models the so-called 
self-shading effect. We give further details of the parameterization of PAR in 
Section 6, and here we only suppose it is a positive function, continuous in P 
for a.e. t, x and measurable in t, x for all P. In order to prove the existence 
result we have to notice that: 

Li is defined on [0, T] x [0, L] x IR, 

< Li{t, X, P) < 1, a.e in [0, T] x [0, L] x IR, 

{t, x) Lj{t, X, P) is measurable, for all P e IR, 

P Li{t,x,P) is continuous, for a.e {t,x) € [0,T] x [0,L]. 

dD 

5. Eventually, let us remark the presence of the advection term u<j-r— in the 

ox 

detritus equations. Detritus, D, sink at a speed of Vd. 



3. Mathematical prelimiriEiries and statement of main result 
3.1. Functional spaces 

In this section we introduce the functional spaces which we use in the remainder 
of this work. All this study is conducted on the open set ]0, L[ and T is a fixed 
time. Throughout this work, concentrations, Cj, are considered as elements of 
the functional space L^(0,L) whose Hilbert space structure is convenient to use. 
However, let us remember that L^(0,L) is continuously imbedded into L^{0,L) 
which is a natural space for concentrations. 
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H and are the separable Hilbert spaces defined by 

H-(L2(o,L))4, 
Hi = (Fi(0,L))4. 

H is equipped with the scalar product 

(C,C) = / y2a{x)Ci{x)dx 

•^0 ^=l 

4 

= ^(Ci, C'i)i2(o,Z,). 

We denote by 1 1 . 1 1 the induced norm on H. 



is equipped with the scalar product 

/.L 4 -.L 4 

(C,C)i = / y2Ci{x)Ci{x)dx + / V 

•^0 i=l -^0 i=l 



dQ{x) dC,{x)^^ 



dx dx 

1=1 i=l 

We denote by ||.||i, the induced norm on H^. 

We will also have to consider the space L°° = {L°°{0, i))**. L°°{0, L) is a Banach 
equipped with the norm 

liailoo =m/{M;|C,(a;)| < M a.e. in (0,L)}. 

Similarly L°° is a Banach space equipped with the norm 

||C||oo= sup ||Ci||cx>- 
j=l,...,4 

Now, if X is a real Banach space equipped with the norm \\-\\x, C{\O^T],X) 
is the space of continuous functions on [0,T] with values in X, equipped with the 
norm, 

||C||c([o,T],x) = sup ||C(i)||x- 

[0,T] 

Similarly -L^(0, T, X) is the space of functions in time with values in X, equipped 
with the norm, 

\\c\\L^,.T.x) = {f\\cm\dt)^'\ 

and L°°(0,T, X) is the space of functions L°° in time with values in X, equipped 
with the norm, 

l|C||i~(o,T,x) = inf{M; \\C{t)\\x < M a.e in {0,T)}. 

C{[0,T],X), L2(o,T,X) and L°°{0,T,X) are Banach spaces. 
We have the useful 
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Lemma 3.1 The imbedding, C L°°, is continuous. 

The imbeddings, c H and c C([0, L], IR"^) are compact. 

Proof. It is a consequence of corollaries IX. 14 and IX. 16 in Brezis ^, and of the 
Rellich-Kondrachoff theorem (Lions and Magenes U . 

H' denotes the dual of H and (H^)' the dual of H^. When H is identified with 
its dual, we have the classical scheme, 

C H = H' c (H^)', 

where each space is dense in the following and the imbeddings are continuous. 
Let us denote by W(H^) the Hilbert space, 

W{H^) = {C e L\0, T, H^); ^ G L^O, T, (H^)}. 

Lemma 3.2 Every C € W^(H^) is a.e equal to a continuous function from [0,T] 
to H. Moreover we have the following continuous imbedding, 

W(Hi)cC([0,T],H). 

Proof. See Dautray and Lions for example □ . 

Moreover, because the injective mapping c H is compact, we know that. 
Lemma 3.3 The identity mapping, W^(H^) c L^(0,T, H), is compact. 

Proof. See Aubin ^ or Lions □ . 



3.2. A preliminary transformation of the system and the bilinear form 

a(i,C,C') 

In order to work with a bilinear form as simple as possible, we start by adding Xd 
to both sides of system (2.1). The value of A > will be fixed in what follows. This 
leads to the equivalent system, for i = 1 to 4: 

= f,{t,x,c) + xCi, ie]o,r], xe]o,L[, ^^^^ 

f (,o)^f (,L) = o, t.M 

Ci{0,x)=C°{x), x€]0,L[. 
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For N, N', P, P', Z, Z' and D, D' e H^O, L), we define 

/■^ ON dN' 

aN{t,N,N') = J d{t,x)^^ + Xj NN', 

dP dP' 

a.it.Z.Z')^l „(,..)-- + A ZZ: 

dD dD' dD 
an{t,D,D') = d{t,x)^^ + v^^D' + \ DD' , 



and 



a(i, C, C) = ajv(i, N, N') + ap{t, P, P') + az{t, Z, Z') + aoit, D, D'). 



Lemma 3.4 For a.e. t e [0,T], a(t, C, C) is a continuous bilinear form on x 
H^. For all C, C e H^, t — )• a(f, C, C) is measurable and there exists a constant 
Ma > such that, 

|a(i,C,C')| < M„||C||i||C'||i, VC,C' e H\ 

For a fixed A, A > 777-, there exists a constant cq > such that, 
2 do 

a(t,C,C) > co||C||f, Vte[0,T], VC e 
Proof. The proof for this is classical and it is omited □ . 



3.3. The reaction terms and the nonlinear operator G 

In this paragraph we show that the reaction terms of the NPZD-model enable us 
to define a continuous operator G on L^(0,T, H). 

Lemma 3.5 The reaction terms fjsr, fp, fz and fo defined in Section 2 have the 
following properties: 

(PI) For a.e. {t,x) e [0,T] x [0,L], and all C e IR^, 

|/jv(i,a;,C)| < - 7) + t)|P| + (m. + t)|Z| + (/xd + r)|£>|, 

\fp{t,x,C)\ < (Mp(l-7)+mp + r)|P|+5,|Z|, 

\fz{t,x,C)\ < {{ap + ad)gz + rriz + Hz + t)\Z\, 

\fD{t,x, C)| < (((1 - Op) + ad)gz + mz)\Z\ + mp\P\ + (/Xd + t)\D\. 



(P2) The function f{t, x, C), defined from [0, T] x [0, L] x IR^ IR*, is measurable 
in {t, x), for all C e IR^, and is continuous in C, for a.e. {t, x) e [0, T] x [0, L] . 
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Proof. The proof is straightforward and uses the comments of Section 2.2 □ . 

We now define a function g{t, x, C) = i{t, x, C) + AC from [0, T] x [0, L] x IR* ^ 
IR* and a nonUnear operator, G, by: 

GC = g{t,x,C{t,x)), {t,x) e [0,T] X [0,L]. 

Proposition 3.1 The operator, G, is well defined from L^{0,T,'H.) to itself. There 
exists a constant Mg > 0, depending only on the parameters of the model, such that, 
for all C e £2(0,T,H) 

l|GC||£,2(o,T,H) < Afg||C||i2(o,T,H)- 

The operator G is continuous on 1/^(0, r,H). 

Proof. Let C e L2(0,T,H) and t e [0,T]. Prom point (PI) of lemma 3.5 we 

obtain, 

\\Gcm' = 

I \fN{t, x, C{t, x)) + XN{t, x)f + \fp{t, x, C{t, x)) + XP{t, 
Jo 

+\fz{t, X, C{t, x)) + XZ{t, + \fn{t, x, C{t, x)) + XD{t, x)\^dx, 

WGcrn' < 

ciei(i|p(t)iii.(„^^) + ii^wiii.(o,^) + iii?wiii.(ox) + \mmUo,L)) 

+CteMPmh(o,L) + \\m\hiO,L) + II^WIIi^(0,L)) 

+cte3i\\Z{t)\\l,^,^^^) 

+cte4(||Z(t)||i,(„^^) + ||P(<)lli.(o,L) + \\m\\Uo,L))^ 
||GC(t)||2<Af|||C(t)|^ 
and integrating on [0,T], 

||GC||z,2(o,T,H) < Afg||C||i:,2(o,T,H)- 

Prom point (P2) of lemma 3.5 we know that the function, 

g(i,a;,C) =f(i,a;,C) + AC, 

from [0, T] x [0, L] x IR"* IR*, satisfies the conditions of Caratheodory and by the- 
orem 2.1 page 22 of Krasnosel'skii we know that the operator G is continuous □ . 



3.4. Variational formulation 

We can now write the definition of a weak solution to system (2.1), 
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Definition 3.1 C e T4^(H^) is a weak solution of system (2.1) if 
V,/.eHi, (^,0)+a(i,C,<^) = (GC,</.), 

in the X>'(]0,T[) sens, 
and C(0) = 0°. 

and state the main result of this paper, 

Theorem 3.1 Let C° € H. There exists a weak solution to system, (2.1). Further- 
more, if N^jP'^jZ^ and D° are positive then N,P,Z and D are positive for a.e. 
t&[0,T]. 

The proof is given in the next two sections. 
4. Existence 

The existence result is obtained in two steps. We first define an approximate prob- 
lem, in which the operator G is approximated by an operator G„. This approximate 
problem is solved using the Schaudcr fixed-point theorem. In the second step we 
let n — >■ 00 to obtain a solution to the initial problem. 

4.1. Step 1: approximated problem 

Let n > be a fixed integer and g„ be defined by, 
g„: [0,T] X [0,L] X IR^ ^ ffi,^ 

^+ n\{s{t,X,C))i\ 

Define the nonlinear operator, G„, by: 

G„C = g„(t, X, C{t, x)), {t, x) e [0, T] X [0, L]. 

Proposition 4.2 The operator, G„, is well defined from 1/^(0, r,H) to itself and 
there exists a constant, Mg > 0, such that for all C e i^(0, T, H), 

||G„C||i:,2(o,T,H) < Afg II C| I L2(o,T,H)- 

The operator G„ is continuous on i^(0,T, H). 

For all C e L^(0,T, H), we also have the estimation, 

l|GnC||i2(o,T,H) <2n\/lT. 

Proof. Let C € Z/^(0,T, H). Prom the definition of g„ and from proposition 3.1 
we obtain, 

l|G„C||i2(g 2-,H) < l|GC||i2(g_T_H) < -^g I I C | | ^2 (Q^-p^H) • 

The estimation | |G„C| |i2(o^T,H) ^ 2n^ LT is also derived easily from the choice 
made to define g„. 
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As in proposition 3.1, g„ satisfies the conditions of Caratheodory and G„ is con- 
tinuous on ^2(0, T, H) □ . 

We now seek a solution to the approximated system and show that such a 
solution is a fixed-point of the operator & defined in the next proposition. 
Proposition 4.3 Let C be a fixed element o/L^(0,T, H) and let C° e H. There 
exists a unique solution to the problem: 
find C e W{U^) such that, 

y^eU\ (^,</.)+a(t,C,</.) = (G„C,<^), 

in the V'{]0,T[) sens, 
and C(0) = CO. 

This solution defines an operator & on L^(0,T, H), @C = C. 

Proof. Since the problem is linear in C and G„C is fixed in L^(0, T, H), the proof 
is classical (e.g. Dautray and Lions, □ . 

To insure that G has a fixed point, we show that the Schauder fixed point 
theorem can be applied. 

Lemma 4.6 The operator 6 is continuous on i^(0,T, H). 

Proof. Let and e L^{0,T,H). and C^, the associated solutions to the 
problem of proposition 4.3, satisfy, 

(^(C^ - C^), </.) + a{t, - <P) = (G„Ci - G„C2, cl>). 

Taking </> = — as a test function, integrating on [0,t], using the coerciveness 
of a and Cauchy-Schwarz inequality, we obtain. 



L 



< f \\GnC\s)-GnC\s)\\\\C\s)-C^{s)\\ds. 
Jo 

As C^(0) = C^(0) = C° we obtain using Young's inequality, 

\\C\t) - C2(t)||2 + f 2co||C^(s) - C^s)\\lds 



f-t 

< 



[ -\\Gr,C\s)-GnC^is)\\^ds + a [ \\C\s) - C^s)\\^ds, 

Jo OL Jo 



and with a = 2co, 



\C\t) - C'{t)f < r -L||G„Ci(s) - G„C2(^)||2^^_ 
Jo 
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Eventually, we obtain, integrating on [0,T], 

' 1 ' 2 12 / T ~ I "2 

||9C — 9C ||l2(o_7-_h) = ||C — C ||l2(o,t,h) < Y 2^11*^" ~ C1„C ||l2(o,t,h) 
and © is continuous as G„ is □ . 

Lemma 4.7 The operator 6 maps L^(0,T, H) in the ball 



,2LTn^ . „^0||2x 



B = {Cei:^(0,T,H),||C|U2(o,T,H) < yT(^— + ||C"||^)}. 

In particular, we have, 0(B) C B. 

Proof. Let C € iy^(0,T, H). C, the solution to the problem of proposition 4.3, 
satisfies, 

{j^C,ct>) + a{t,C,^) = iGnC,(t>). 

Taking </> = C as a test function, integrating [0,t], using the coerciveness of a and 
the Cauchy-Schwarz inequality we obtain, 

Jo ^1"^^'^"' + co\\Cis)\\lds < 1^ ||G„C(s)||||C(s)||ds, 
with Young's inequality, 

\\C{t)\f+ f 2co\\C{s)\\lds < f -\\GnC{s)\\^ds + a f \\C{s)\\^ds + \\CY, 
Jo Jo Oi Jo 

with a = 2co and as / ||G„C(s)|p(is < 4Lrn^, we have 
Jo 

llcw|p<^ + llc°||^ 

integrating once more on [0, T] we obtain. 



Co 

□ . 

Lemma 4.8 The operator 6 is compact. 

Proof. Let i? be a bounded set in L^{0, T, H). Let us show that Q{B) is bounded 
in VK(H^). Let C e B C L^(0,T, H) and let C be the associated solution to the 
problem of proposition 4.3. As in the proof of lemma 4.7 we obtain 

\\C{t)\f+ f 2co\\C{s)\\lds < f -\\Gr,C{s)\fds + a f \\C{s)\\^ ds + \\CY ■ 

Jo Jo 01 Jq 
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Taking a = cq this time, we obtain 



I 



Co Co 

and eC is bounded in L'^{0,T,ll^). 
Moreover we obtain 

V<^eH\ {^,c^)+a{t,C,ct>) = {GnC,^). 

From lemma 3.4, 
r/C 

|(— ,.^)| < M„||C||i||</.||i + ||G„C||||<^|| < M„||C||i||^||i + 2r^^/I||.^||l, 

r W^WUryds < 2 r{4Ln^ + M^\\C\\l)ds. 
Jo dt Jq 



and 



And therefore, ||^^||l2(o (hi)') is bounded in L^(0,T, (H-^)'). 

The range of 6 is in VK(Hi), from lemma 3.3, the injection W{\1^) C L'^{Q,T,'H.) 

is compact, and this concludes the proof □ . 

It is now possible to state the main result of this section, concerning the existence 
of weak solutions to the approximated problem. 

Theorem 4.2 Let n > be a fixed integer. Let C'' G H. There exists a solution, 

C„, to the problem: 

find C e W{lV-) such that, 

V</.eHi, (^,</>)+a(t,C,</>) = (G„C,0), 

in the V'{]0,T[) sens, 
and C{0) = C°. 

Proof. From lemma 4.6, 4.7, 4.8 and the Schauder fixed-point theorem, the oper- 
ator Q has a fixed point, which is the solution sought □ . 



4.2. Step 2: letting n ^ oo 

We now pass to the limit as n ^ oo, in the equations, 

dC 

V</)eHi, {^{t),cl>)+a{t,Cn,(t>) = {GnCnit),cl>), ^4 4^ 

C„(0) = CO. 

This is achieved in two steps: 



6 Well-posedness of a coupled biological-physical model for the upper ocean 

a) a priori estimations on the sequence C„, 

b) extraction of subsequences and letting n — )• oo. 
a) estimations. Let us show that: 

(a.l) the sequence (C„)„>o is bounded in L°°(0,T, H), 
(a.2) the sequence (Cn)n>o is bounded in L'^{0,T,H^), 

(a.3) the sequence (~^~)">o is bounded in L^{0,T, (H^)'). 
Taking, C„ as a test function in (4.4) we obtain, 

2 ^l|Cn||^ + a{t, C„, C„) = (G„C„, C„), 



or, 

1 d 



C„|p + co||C„||f < ||G„C„||||C„|| < Mg||C„| 



2dt 

and integrating on [0, t], we obtain 



||C„||2 + 2co/ ||C„||frfs<2Mg / ||C„||2ds+||C' 
Jo Jo 

Equation (4.5) gives 



t 

2^o)i/r / iir~" i|2j„ I iir~<0||2 



\\Cnr<2Mg / ||C„||^ds+||C 
Jo 

Using Gronwall lemma, we have 

||C„(i)ll'< l|C°||'exp(2M,T), 
and the sequence (C„)„>o is bounded in L°°(0,T, H). 

Equation (4.5) also gives 

f\\Cn\\lds<^ f\\Cn\\^ds + ^\\C 
Jo Co Jo 2co 

and with (4.6) 

r ||C„||?ds < ^||C0||2exp(2M,T) + -L||C 
Jo Co ZCo 

Therefore the sequence (C„)„>o is bounded in L^(0,T, H^). 



0||2 
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dC 

Let us now give an estimation for the sequence, ( ^^" )">o- have 



" dt 
and therefore 

..dC 



dC 

",.^)| < |a(t,C„,.^)| + |(G„C„,< 



,</>)! <M„||C„||i||<^||i+M3||C„||i||<^||i, 



" dt 
that is to say 

£\\^\\h^y<2{M'^+M^) £\\Cn\\l 



M T 1 
< 2(M2 + M2)(^||C°|Pexp(2M„T) + — ||C°ii2^ 
Co 2co 



and the sequence {—^)n>o is bounded in L^{0,T, (H^)'). 



b) passing to the limit. Let us first recall that V{]0, T[, H^) C W{U^). There- 
fore e and \/ip e V{]0,T[), we have ij; = (fxg) (p e L'^{0,T,li^) and 

''''' GL\0,T,{1i^y). 



dt 

dC 

b.l) The term from (a.3), we are able to extract from the sequence 

dC 

( ^^ )n>o a subsequence (denoted in the same way) converging to some 
h in i^(0, T, (H^)') weak star, that is to say, for all ^ G and all 

vpe2?(]o,r[). 



n->-oo Jq dt Jq 



Moreover, by definition, we obtain 



Prom (a. 2), we are able to extract from the sequence (C„)„>o a subse- 
quence (denoted in the same way) converging to some C in L^{0,T,'H.^) 
weak. Therefore, for all (peH^ and all ip G V{]0,T[), 

lim - r {Cr.,cl,)^ds = - r {C,4>)'^ds, 
Jo Jo dt 

dO 

andh = ^ in 1.2(0, T, (Ri)'). 
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b.2) The term a(t, C„,0): from (b.l), we can suppose that the sequence 
(C„)„>o converges to C in L^(0,r, H^) weak. Therefore the sequence 
{dxCn)n>o converges to dxC in L^(0,T, H) weak. Then, for all 4> £ 
and all ip e V{]0,T[), 

lim / a{s,Cn,(l))ipds = / a{s,C,(j))ipds. 
Jo Jo 

b.3) The term (G„C„,(/)): from (a. 2), (a. 3), and from the compacity of 
the injection W^(H^) — )• L^(0, T, H), we can suppose that the sequence 
(C„)„>o converges to C in L^(0,T, H) strong. Therefore, each C„ i, 
i = 1, ...4, converges to Ci in L^(0, T, L^(0, L)) strong. From the inverse 
Lebesgue theorem (Brezis, ^, theorem IV. 9. page 58), we can suppose 
that: 

(b.3.1) the sequences (Cn,i)n>o, i = 1, •••4, converge to Ci a.e. in]0,T[x]0, L[. 
(b.3. 2) for i = 1 to 4, |C„,i| < hi, Vn > 0, a.e. in ]0,r[x]0,L[ and hi e 
L2(0,r,L2(0,L)). ' 

As gn,i(i, X, C) is continuous in its third variable, we deduce from (b.3.1) 
that \/4>i e H'^{Q,L) and ^ip e ©(]0,T[), 

Un,i{t,x) = gn,i(t,X, Cn(t,x))4>i{x)(p{t) 
_^ gi{t,X,C{t,x))^i{x)<fi{t), 

n — > OO 
a.e. in ]0,r[x]0,L[. 

Moreover, from lemma 3.3 and with (b.3. 2) we have, 

4 

\un,i\ < MiiJ^hiMiM e L\]0,T[x]0,L[), 
1=1 

where the Mi arc constants. Thus, from the Lebesgue theorem on dom- 
inated convergence, we obtain, 

pT pL pT pL 

lim / / Un i{t,x)dxdt = / / gi{t,x,C(t,x))(l)i(x)ip(t)dxdt, 

"-*°°Jo Jo ' Jo Jo 

and finally, for all <j)eU^ and all (peV{]0,T[), 

lim / (G„C„,0)^= / {GC,(P)^. 
"^°°io Jo 

This concludes the proof of the existence of weak solutions to the one-dimensional 
NPZD-model. 



5. Positivity 
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In this section we prove the second part of Theorem 3.1: if initial conditions 
N^, P°, and are positive then solutions to the one-dimensional NPZD-model 
are positive for a.e. t G [0,T]. To prove this, we need to treat each of the four 
equations seperately, in detail, and in a convenient order. We first show that Z and 
P are positive. Next we show that D is positive using the fact that Z and P are 
positive. Finally, as Z, P and D are positive we obtain the positivity of N. 

• Let us recall that for aU C G and all t € [0,r], aN{t,N,N) > 0, 
ap{t,P,P) > 0, az{t,Z,Z) > and az{t,Z,Z) > 0. 

• Z is positive: 

Let C be a weak solution to the NPZD-model. Let us take 

—Z~ = —max{0,—Z), 

as a test function. Since, 



and 

az{t, Z{t), -Z-{t)) = az{t, Z (t)- , Z {t)- ) , 

we obtain 

^|lI^(i)-||i^(o,L) + = -(5z(C(i)),^(i)-)- 

Let us detail the term (gziC), Z^). 

{gz{C),Z-)L2^o,L)= {ap^^ZZ-+ad^^^ZZ--m,ZZ- 
-H,ZZ-)t^o,i] + {-tZZ-)\i^l[ + ^ZZ-. 
As ZZ- = -{Z-f, we have 

\2 HzJ^ x/r,_\2 



)iz-Y - {a,-^^){z-y + m^z-y 



and 

32 



igz{C),Z-)L2^o,L)> li-i^Pj^^KZ- 
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or 

Thus ^ 

-{gz{C),Z-)L2(^o,L) < I {g,{ap + ad) + X){Z-f), 

Jo 

As az{t, Z~,Z~) > 0, we obtain 

^ll-^"lli=(o,L) < Hgziap + ad) + A)||Z-|||2(o,L)- 
Integrating this inequahty on [0, t] and using Gronwall's lemma, we obtain 
< ll^"(0)||i2(o,i) exp(2(ff,(ap + a^) + X)t). 

Therefore Z is positive. 
P is positive: 

In the same manner, let us examine the term (5p(C), P~)i2(Q^). 

(5P(C), P-)lho,l) = j\np{l - ^)LiLmPP- - ijT^ZP-) 
-nipPP-) l]o,i] + {-tPP-)1]i^h + XPP-, 

= j\-l^p{-^ - l)LjL^{P-f - ij^^ZP-) 
+mp(P-)2) l]o,i] + (r(P-)2) - A(P-)2, 

> / {-Hpil-j)LjLN{P-f 
Jo 

and 

-(ffp(C),P-)i2(o,i) < / f,p{l-^)LiLN{P-f 

Jo 

— X 1 

The function x i-)- j bounded by — = on IR. 

Z/jv and Lj are bounded by 1. 
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Z{t), is a solution to the NPZD-model and therefore belongs to H^{0,L) C 
L°°(0,L). Hence we have, Vt, Z{t) < \\Z{t)\\oo, and 

-{gp{C),P-U2^o,L) < (A + Mp(l-7)+5z^||^(i)lloo)||P-Wlli2(o,L)- 
We conclude in the same way to obtain 

\\P~it)\\h(o,L) < ll^"(0)lli^(o,L)exp(^ 2(A + Mp(l-7) 

+9z^J\Z{s)\\oo)ds). 

D is positive: 

(gniC), D-U.^o,L) = - a,){^^ZD-) - a,{^^ZD-) 

+mpPD- + m,ZD- + MD-f) ^]o,i] 
+ iT{D-f)l]i,L[ + XDD-. 
Because P, Z and D~ are positive, we obtain 

{goiC), D-)r.^o,L) > £ V^] - ^(^")'. 

-{gniC), D-)i.(o,L) < j\-aa{^^Z{D-f)) l^^A + A(I)-)^ 

-(5z5(C),D-)^.(o,^) < j\aa[g,^\\Z\\^){D-f)l^^^i^+\{D-)\ 

-(5j3(C),L'-)^2(o,i) < (A + ad(5z^^||^||oo))||-D"||i2(o,r). 
Hence 

|ll^-|li^(o,L) < 2{X + aa{g.^\\Z\U)\\D-\\l.^^o,LY 
and we can conclude. 
A'' is positive: 



{9 



N 



(C),7V-)i2(o,i) = / {-f,p{l-^)LiLNPN- +f,,ZN- 
Jo 



+fidDN-) Ijo,;] + (t(P + Z + D)N-) + AA^iV". 
Because P, Z, D and A''" are positive, we have 

(ffjv(C), iV-)i2(o,L) > ^\-^ip{l - ^)LjLr,PN-)l^o,i] - A(A^-)^ 
Jo 

-{gN{C), N-)l^o,l) < 1^ - l)Li^-^^P{N-f) l^o,i] + X{N-f. 
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Once again we can conclude and the proof of theorem 3.1 is complete. 

Hence, if initial concentrations arc positive then concentrations are always pos- 
itive and both models, with or without absolute values in the nonlinear terms, are 
equivalent. 

6. Existence, positivity and uniqueness for different Gp, Gdi Li and 
zooplankton mortality formulations 

Functions used to parameterize biological fluxes such as zooplankton grazing on 
phytoplankton, Gp, or on detritus, Gd, light limited growth rate, Li or zooplank- 
ton mortality (which is a constant, m^, in our model), vary from one modeling study 
to another. One can wonder if the existence result still applies with these different 
formulations. To answer this question it should be noticed that tlie key argument 
used in the proof is the fact that the nonlinear reaction terms allow us to define 
a nonlinear continuous operator G satisfying ||GC||^2(g_2i jj) ^ -^3I|C||l2(o,t,h)- 
Therefore, as all the functions listed in Table 2, found in the literature, are contin- 
uous and bounded on ffi,"*" or (IR"'')^, the existence result stays correct. Positivity 
can also easily be checked for all these different formulations. It should however be 
mentioned that some studies use a quadratic zooplankton mortality term which can 
not be treated with the method we propose. 

Let us now concentrate on the question of the uniqueness of weak solutions to 
the one-dimensional NPZD-model. 

In order to prove uniqueness we need the nonlinear reaction terms to satisfy a local 
Lipschitz condition which was not needed to obtain the existence result. To ver- 
ify that such a condition holds we examine in some details the optical model from 
which the PAB.{t, .x, P) term and consequently the Li{t, x, P) term are calculated. 
In the different equations Li{t, x, P) allways appears in the product form PLi{t, x, P) 
Concerning this product the desired local Lipschitz condition reads as follows: 

for all {t,x) G [0,r] X [0,L], and all F, P € [0,-hoo[, 

\PLi{t,x,P)- PLi{t,x,P)\<Ki{P,P)\P - PI (6.7) 

where Ki is a continuous nonnegative real-valued function which is increasing in 
each variable. 



In the optical model we considered, two different wavelengths are taken into account 
and the absorption coefficients depend on the local phytoplankton concentrations: 

PAR{t,x,P) = Q(i)(exp(-(fcgo + kgp{^^^)^^)x) 
+ exp{-{kro + krp{ )''')a;)). 

Q(t) is proportional to the irradiance intensity hitting the sea surface at time t. 
Parameters are given in Table 3. Let us suppose that Q{t) G L°°{0,T) and that 
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Table 2: Different parameterizations found in the literature. All parameters are 
positive constants. 



Z grazing on 
P, Gp 




this study and e.g. 
Fennel et al. 


k, + P^ + £>2 


e.g. Losa et al. 


9zrP^ 


e.g. Fasham et al. 

10 


k^{rP + (1 - r)D) + rP^ + (1 - r)D'^ 


Z grazing on 
D, Gd 




this studv and e Ef 
Fennel ei a/. 




T f 7 21 

e.g. Losa al. 




e.g. Fasham aZ. 

10 


k^{rP + (1 - r)D) + rP2 + (i _ r)D^ 


light limited 
growth rate, 
Li 


1 - exp(-P^P(i, X, P)/kpar) 


this study and e.g. 
Levy et al. 


VpaPAR{t,x,P) 
{v^ + a'^PAR{t,x,Py)y^ 


e.g. Spitz et al. 


Z mortality 




this study and e.g. 
Levy et al. ^"^ 


rUzZ 
k + Z 


e.g. Losa et al. 
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Q{t) > 0. Even though the exponents Ig and Ir satisfy < Ig < 1 and < < 1, 
an easy calculation of the derivative, --j^{PLi{t,x,P)), shows that with such an 
optical model property (6.7) is satisfied with: 

Ki{P,P) = 1+Mk^:(fc / (i^)'»(max(P,P))'«+fc,A(^)''^(max(P,P))''-). 

f^par ipgic 'pg'c 

In the literature PAR{t, x, P) is often parameterized by, 

PAR{t, X, P) = Q{t) exp(-(fci + k2P)x), 

whore &i and ^2 are positive constants. With this simpler formulation property 
(6.7) is clearly satisfied. 

The following two lemmas give the local Lipschitz property satisfied by all four 
reaction terms of the NPZD-model. 

Lemma 6.9 The nonlinear reaction terms gjq, gp, gz and go satisfy: 
for all {t,x) € [0,r] X [0,L], and all C,C € (IR+)^ 

\gN(t,x, C) - gN{t,x, C)| < Kn{P. P){\N - N\ + \P - P\ + \Z - Z\ + \D - D|), 
\gp{t, X, C) - gp{t, x,C)\< Kp{P, P){\N - N\ + \P - P\ + \Z - Z\ + \D - D\), 
\gzit, X, C) - gz{t, x, C)| < Kz{Z, Z){\N - N\ + \P - P\ + \Z - Z\ + \D - D\), 
\gn{t, x, C) - goit, x, C)| < Kn{Z, Z){\N - N\ + \P - p\ + \Z - Z\ + \D - D\), 

where Kn, Kp, Kz, Kd are continuous nonnegative real-valued functions which 

are increasing in each variable. 

9 

X X 

Proof. Functions l{x) = and g{x) = ^ are continuously differentiable 

on [0,+cxd[, and 

\l'{x)\ < -1, \g'{x)\ < 



S's/k;^ 

Therefore I and g are Lipschitz continuous. 

It is clear that 

\gN{t, X, C) - gN{t, x,C)\ < - ^)\Lj{t, x, P)PLn{N) ~ Lj{t, x, P)PLn{N)\ 

+{fi, + t)\Z -Z\ + [lid + t)\D- D\ 
+t\P - P\+X\N - N\. 

Now since 

\Li{t,x,P)PLN{N) - Li{t,x,P)pLN{N)\ = \Li{t,x, P)P{Ln{N) - Ln{N)) 

+LNiN){Li{t, X, P)P) - Li{t, X, P)P\, 

we have 

\Li{t, X, P)PLN{N)-Li{t, X, P)PLn{N)\ < max (P, P)\N-N\+Ki{P, P)\P-P\. 
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We then define 

Kn{P,P) — max (A + — 7)- — max (P, P), r + Hp{l — j)Kj(P, P), fx^ + (Xd + t). 
Kp, Kz and Kd are obtained in the same way □ . 

Lemma 6.10 Fort G [0,T] and for positive C{t),C{t) G H^, there exists a con- 
stant Lao, depending on \\C{t)\\ 00 a'nd\\C{t)\\ao, such that the operator, G, satisfies, 

||GC(i) - GC(i)|| < L^\\C{t) - C{t)\\. 

Proof. Prom lemma 3.1, c L°°. Prom lemma 6.9 we have, 

||GC(t)-GC(t)||2 

= / \gN{t, X, C{t, x)) - gN{t, X, C{t, x))\'^ + \gp{t, x, C{t, x)) - gp{t, x, C{t, x))\'^ 
Jo 

+ \gz{t, X, C{t, x)) - gz{t, x, C{t, x))\'^ + \gD{t, x, C{t, x)) - goit, x, C{t, x))\'^dx, 

< cte{Kr,{\\P{t)\\^, \\Pit)\Uf [ \\N{t) - NmUo,L) 

+iip(i) - p(t)iii.(o,^) + \\zit) - m\\Uo,L) + mt) - m\\Uo,L) ] 
+... 

< ioo(||C(t)||oo, ||C(i)||oo)'||C(i) - C(t)||2. 
□ . 

Elementary calculations show that functions of Table 2 are continuously dif- 

ferentiablc on IR"*" or (IR^)^, with bounded first derivatives. Therefore they arc 
Lipschitz continuous and the uniqueness result presented below also holds for these 
formulations. 

Proposition 6.4 The weak solution to the one- dimensional NPZD-model, C e 
W(H^), is unique. 

Proof. Let us suppose that there are two solutions and e W^(H^). They 
satisfy 

V<^ e U\ (^(C^ - C^), cp) + a{t, - C\ 4>) = (GC^ - GC^ </.), 

Let us choose (/) = — as a test function. Using the coerciveness of a and the 
Cauchy-Schwarz inequality, we obtain 

i|||Ci(i)-C2(i)ll' + co||Ci(t)-C2(t)||?<||GCi(i)-GC^(t)||||Ci(i)-C2(t)||. 
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With Young's inequality, we obtain 

ll\\C\t)-C'{t)\\' + co\\CHt)-C'{t)\\l 



with a = 2co, 



|||Ci(t) - C\t)\\' < ^J\GC\t) - GC'{t)\\\ 



From lemma (6.10) we obtain 

j^\\C\t) - C\t)f < ±LUmc\t) - C\t)f. 
Thus, integrating on [0, t] and using Gronwall's lemma 

||C^(i) - C\t)\f < \\C\0) - C^(0)|r exp(^* ^^-(^)'^^)- 
This concludes the proof □ . 



Table 3: Optical model parameters 



parameter 


name 


value 


unit 


Rcdficld ratio C:N 


rd 


6.625 




contribution of Chi to absorbing 


'''pg 


0.7 




pigments 








carbonc:chlorophyll ratio 




55 


mgC.mgChla~^ 


water absorption in red 


h-pQ 


0.225 




water absorption in green 


kgo 


0.0232 




pigment absorption in red 


krp 


0.037 


.{mgChl.m~^)~^^ 
m~^.{mgChl.m~^)~'<' 


pigment absorption in green 


kgp 


0.074 


power law for absorption in red 




0.629 




power law for absorption in green 




0.674 





7. Conclusion 

We have presented a qualitative analysis of a one-dimensional biological NPZD- 

modcl. This model describes the evolution over time and space of four biological 
variables, phytoplankton, zooplankton, nutrients and detritus. The only physical 
process which is taken into account is vertical diffusion and the biological model 
is imbedded in a physical turbulence model which we did not give explicitly but 
appeared as a space and time-dependent mixing coefficient. The model's equation 
for detritus also contains an advection term which represents the sinking of detri- 
tus with a constant speed. All four variables interact through nonlinear reaction 
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terms which depend on space and time through the action of light and present a 
discontinuity in the space variable at a particular depth. 

We have formulated an initial-boundary value problem and proved existence of 
a unique weak solution to it. Furthermore, a detailed investigation of the reaction 
terms enabled \is to prove positivity of the solution. This is biologically important 
since variables represent concentrations which should always be positive quanti- 
ties. We have also shown that the result still holds if different parameterizations of 
biological processes found in the biogeochemical modeling literature are used. 

The analysis conducted in this paper is a necessary first step towards the inves- 
tigation of qualitative properties other than positivity which might be of interest. 
For example Boushaba et al. ^ deal with the problem of determining the asymp- 
totic behavior of solutions to their phytoplankton model. One could also wish to 
investigate the bifurcational structure of the NPZD-model even though the com- 
plexity of the analytical formulation of the equations might constitute a difficulty. 
This type of study could help the understanding of the modifications of evolution 
of the NPZD system under minor changes in the values of parameters reported in 
Edwards ^ and Faugeras et al. 

Eventually we would like to point out that the analysis we presented can easily 
be extended to models containing any mimbcr, n, of biological variables as long as 
the nonlinearities allow us to define a continuous nonlinear operator G satisfying 
||GC||l2(o_t,(l2(o,l))") < Ai"g||C||L2(o,T,(L2(o,L))")- Howevcr, in such more complex 
models, the question of positivity seems to be delicate as equations have to be 
treated one after the other, and the right order has to be found as the positivity of 
some variables can depend on the positivity of others. Let us also mention that the 
analysis can be extended to three-dimensional models in which not only mixing coef- 
ficients but also velocities, calculated by an ocean circulation model, are included in 
the system of partial differential equations constituing the biological model. These 
velocities can be included in the bilinear form, a{t, ., .), as Vd, the detritus sedimen- 
tation speed, is in the formulation of the initial-boundary value problem we studied. 
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